Derivation of the Lattice Boltzmann Model for Relativistic Hydrodynamics 
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A detailed derivation of the Lattice Boltzmann (LB) scheme for relativistic fluids recently proposed 
in Ref. [lj, is presented. The method is numerically validated and applied to the case of two quite 
different relativistic fluid dynamic problems, namely shock-wave propagation in quark-gluon plasmas 
and the impact of a supernova blast-wave on massive interstellar clouds. Close to second order 
convergence with the grid resolution, as well as linear dependence of computational time on the 
number of grid points and time-steps, are reported. 
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I. INTRODUCTION 



Relativistic fluid dynamics plays a major role in many 
fields of modern physics, e.g. astrophysics, nuclear and 
high-energy physics and, lately, also in condensed mat- 
ter. The dynamics of such systems requires solving 
highly nonlinear equations, rendering the analytic treat- 
ment of practical problems extremely difficult. There- 
fore, several numerical methods have been developed, 
based on macroscopic continuum description and 
kinetic theory Very recently, a new Lattice Boltzmann 
(LB) scheme for relativistic fluids has been proposed, and 
numerically validated for two rather different relativistic 
applications, shock waves in quark-gluon plasmas and 
blast waves from supernova explosions impinging against 
dense interstellar clouds [1]. This fills a missing entry 
in the remarkably broad spectrum of LB applications 
across most areas of fluid-dynamics, including quantum 
fluids @ . While a quantitative assessment of its practi- 
cal impact on relativistic fluid dynamics must necessarily 
await for a long and thorough validation activity, work in 
Ref. [l| provides robust indications that the relativistic 
LB (RLB) stands concrete chances of carrying the recog- 
nized advantages of LB schemes for classical fluids, over 
to the relativistic context. We refer primarily to mathe- 
matical simplicity/computational efficiency, especially on 
parallel computers [7J, and easy handling of complex ge- 
ometries. 

In this paper, we present an extended version of our 
previous work [lj. First, we provide full details of the 
analytical and numerical formulation leading to the rela- 
tivistic LB scheme, including the asymptotic (Chapman- 



Enskog) analysis of the continuum fluid-dynamic limit, 
starting from the kinetic level. The numerical validation 
for the case of quark-gluon plasmas is also extended in 
such a way as to probe the convergence/accuracy of the 
RLB schemes as a function of grid resolution. Like for 
classical fluids, second-order convergence and linear scal- 
ing of CPU time with number of grid points and time- 
steps is found. Moreover, the application to supernova 
blast waves is explored in more detail, by investigating 
the effect of increasing Lorentz factors on the space-time 
distribution of the density and pressure fields. Here, the 
numerical simulation of the relativistic flow past a dense 
inter-stellar medium (massive clouds) provides a clear in- 
dication that sweeping of interstellar matter across the 
cloud becomes appreciable only for relativistic beta fac- 
tors (3 > 0.5. 



II. THE BASIC IDEA 

The procedure developed [l[ was prompted by two sim- 
ple observations, i) the kinetic formalism is naturally co- 
variant/hyperbolic, ii) being based on a finite- velocity, 
discrete (beam) representation of the kinetic distribution 
function, standard lattice Boltzmann methods naturally 
feature relativistic-like equations of state, whereby the 
sound speed, c s is a sizeable fraction of the speed of 
light c, i.e. the maximum velocity of mass transport 
(c s /c = K, with 0.1 < K < 1). Based on the above, 
and choosing the lattice speed close to the value of the 
actual light-speed (at variance with standard LB appli- 
cations) ci=5x/6t ~ c, the LB mathematical framework 
allows the relativistic extension developed in our previous 
work. The standard LB reads as follows: 



fi(x + m\ t + St)- fi(x; t) = -uSt {fi ~ ft q ) 



(1) 
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where fi(x; t) denotes the probability of finding a particle 
at lattice site x and time t, moving along the direction 
pointed by the discrete velocity Cj. The left-hand side 
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is readily recognized as an exact lattice transcription of 
the free-streaming term (d t + u V a )/ of the continuum 
Boltzmann equation, where Latin index a labels the spa- 
tial coordinates and repeated indices are summed upon. 
Being naturally covariant, this term goes virtually un- 
changed over to the relativistic context. 

The right-hand-side, on the other hand, is a discrete 
version of the collision operator, here taking the form of 
a simple relaxation around a local equilibrium f^ q , on 
a timescale r = 1/co. The local equilibrium encodes 
the symmetries/conservation laws governing the ideal 
(non-dissipative) fluid regime, namely mass-momentum- 
energy conservation and Galilean invariance. While 
molecular details of the collisional processes can be safely 
foregone, these conservation properties must necessarily 
be preserved in the lattice formulation. 

The discrete local equilibrium is usually expressed as a 
local Maxwellian, expanded to second order in the local 
Mach number Ma = u/c s , u being the local flow speed. 
For the case of athermal flows, this takes the form 



(2) 



where (particle mass is taken to unity for simplicity): 
P = Yl & ' 



pu a 



5> 



are the fluid density and mass current density, respec- 
tively. In the above, Wi is a set of weights, obey- 
ing the sum-rules J2i w i = 1 an< ^ J2i w i c ia = c s< an d 
Qiab — CiaCib — c 2 s 5 a b is the projector along the i-th spatial 
direction. It is readily checked that the local equilibria 
fulfill the following conservation rules 

i i 

£ fV c ™ = £ f iCia = pUa ' (3) 

i i 

^ ^CiaCib = p(u a U b + C 2 s S ab ) . 

i 

The first two are the usual mass-momentum conserva- 
tions laws, whereas the latter ensures the isotropy of the 
equilibrium momentum-flux. The latter is crucial to se- 
cure the proper non-linear structure of the Navier-Stokes 
equations, and indeed only specific classes of discrete lat- 
tices fulfill the aforementioned conservation constraints. 
As previously noted, lattice equilibria can be obtained 
by local expansion of the continuum expression of lo- 
cal Maxwell equilibria. In a more empirical way, they 
could also be obtained by matching the local equilibria 
in parametric form, Ae~ BciaUa , to the conservation rules 
Q, thereby fixing the Lagrangian parameters A and B 
in terms of the conserved hydrodynamic fields p and u a . 



The possibility of fixing local equilibria by simply ex- 
panding the local continuum Maxwellian, which is more 
elegant than empirical matching is by no means evi- 
dent. 

In fact, it is strictly related to the well-known property 
of the local Maxwellian to serve as the generating func- 
tion of Hermitc's polynomials H n (here v — v/c s and 
u = u/cg); 



oo 

E 

n=0 



H n (v)u n 



(4) 



Note that the Galilean invariance manifestly encoded 
at the right-hand side through the dependence on the 
magnitude of the relative speed (v a — u a ), can only be 
preserved by including all terms in the Mach-number ex- 
pansion at the right hand side. It is quite fortunate that 
the Navier-Stokes equations only involve quadratic non- 
linearities in the flow field, because this allows to develop 
a consistent lattice hydrodynamic theory by retaining 
only second order terms in the Mach-number expansion. 
A similar line of thinking can also be applied to the rela- 
tivistic equations, with due changes in the mathematical- 
physical details, to be discussed shortly. 

On the other hand, we are not aware of any relativistic 
analogue of the relation Q for relativistic local equilib- 
ria (Jiittner distribution). Because of this, the relativistic 
LB scheme has been devised according to the moment- 
matching procedure discussed above. That is, the local 
kinetic equilibria are expressed as parametric polynomi- 
als of the relativistic fluid velocity j3 = u/c, with the La- 
grangian parameters fixed by the condition of matching 
the analytic expression of the relevant relativistic mo- 
ments, namely the number density, energy density and 
energy-momentum. As anticipated, the possibility of a 
successful matching stems directly from the fact that, 
even in standard (non-relativistic) LB fluids, the sound 
speed c s is of the same order of the speed of light, typ- 
ically c s — c/V3, which is exactly the equation of state 
of ideal relativistic fluids. As a result, |/3| = Ma/V3, 
so that |/3 1 is of the same order as the Mach number 
Ma = \u\/c s . Thanks to this simple, and yet basic prop- 
erty, it is possible to tackle weakly relativistic problems 
in close analogy with the LB theory of classical low-Mach 
fluids, the algebraic details being of course quite different 
in the two cases. 

This permits to carry most of the LB formalism over 
to the context of weakly relativistic fluids, such as quark- 
gluon plasmas generate d by recent experiments on heavy- 
ions and hadron jets [lOl-tlq. as well as astrophysical 
flows, such as interstellar gas and supernova remnants 
[I2H23. 

The RLB scheme is verified through quantitative com- 
parison with recent one dimensional hydrodynamic sim- 
ulations of relativistic shock wave propagation in viscous 
quark-gluon plasmas 21], and also applied to the three 



dimensional case of a blast-wave, produced by a super- 
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nova explosion, colliding against interstellar massive mat- 
ter, e.g. molecular gas |17| . 

Being based on a second-order moment-matching pro- 
cedure, rather than on a high-order systematic expan- 
sion in /3 of the local relativistic equilibrium (Jiittner) 
distribution [5j, the RLB is limited to weakly relativis- 
tic problems, with ~ 0.1. Note in fact that, un- 
like the continuum Maxwellian, polynomial expansions 
are positive-definite only for Mach-number (relativis- 
tic (3) below a given threshold, typically Ma ~ 0.3. 
However, by introducing artificial faster-than-light par- 
ticles (numerical "tachyons"), the RLB scheme can be 
taken up to |/3| ~ 0.6, corresponding to Lorentz's fac- 
tors 7 = , 1 = ~ 1.4 fijl. Although still far from 

\/l-|/3| 2 MJ 

strongly relativistic regimes, with 7 3> 10 and higher, 
this Lorentz factor is nevertheless relevant to a host of 
important relativistic fluid problems at wildly disparate 
scales, such as quark-gluon plasmas and relativistic out- 
flows in supernova explosions and possibly even Dirac 
fluids in graphene [22l |. 



III. MODEL DESCRIPTION 

We begin our model description by considering the 
relativistic fluid equations associated with the conser- 
vation of number of particles and momentum-energy. 
The energy-momentum tensor reads as follows [23l [24j j : 
T^ v = P^ v + (e + P)u»u v + tt^, e being the energy 
density, P the hydrostatic pressure and n^" the dissipa- 
tive component of the stress-energy tensor, to be specified 
later. The velocity 4-vector is defined by u M = (7,7/?)^, 
where f3 = u/c is the velocity of the fluid in units of 



the speed of light and 7=- 



The tensor 77^" de- 



notes the Minkowski metric. Additionally, we define the 
particle 4-flow, — nj(l, with n the number of 
particles per volume. Applying the conservation rule to 
energy and momentum, d^T^ = 0, and to the 4-flow, 
d^N^ = 0, we obtain the hydrodynamic equations, 



d t (( C + P) 7 2 - P) + d a (( C + Ph 2 U a ) 



aO 







(5a) 



d t ((e + Ph 2 u b ) + d b P + d a ((e + Ph 2 u a u b ) 

n , . (5b) 

+ dtn ob + d a n ab = , 

for the energy momentum conservation, and 

d t {n-f) + d a (wyua) = , (6) 

for the conservation of particle number. Note that, un- 
like the case of non-relativistic fluids, we have two scalar 
equations, one for the particle number and one for the en- 
ergy (in classical hydrodynamics energy appears as the 
trace of a second-order moment, namely the momentum- 
flux) . To complete the set of equations, we need to 
define a state equation relating at least two of the three 
quantities: n, P and e. 



A. Relativistic Boltzmann Equation 

The above hydrodynamic equations can be derived as a 
macroscopic limit of the relativistic Boltzmann equation. 
For the case of a single non-degenerated gas, and in the 
absence of external forces, this reads as follows [23| : 



J P*o 



(7) 



where = y E ^ , pj is the particle 4-momentum with 

E(p) the relativistic energy as function of the momen- 
tum magnitude p=|^|, E(p) — (p 2 c 2 + m 2 c 4 ) 1//2 . In the 
above, f*=f(x,p*,t) and f=f(x,p,t) denote the distri- 
bution functions before the collision, while /*=/(ic,p!|,, t) 
and f'=f(x,p',t) are the resulting ones after the colli- 
sion. The base of the so-called collision cylinder is de- 
scribed by adtt, with a the differential cross section, O 
is the solid angle, and 



c V c z 



(8) 



is the Lorentz invariant flux 23] , with v and v t the veloc- 
ity of the particles with momentum p and , respectively. 
The right-hand-side of Eq. ([7]) is the collision term, whose 
details fix the value of the transport coefficients in the 
macroscopic equations. Although the collision integral 
can be expressed in terms of the second kind modified 
Bessel functions and numerical integrations [23| , simpler 
expressions have been proposed, along the lines of the 
BGK (Bhatnagar-Gross-Krook) approximation for non- 
relativistic fluids. The first relativistic BGK (RBGK), as 
proposed by Marie [25] . reads as follows: 



W7) = — (f eq -f) 



TM 



(9) 



where f eq is a local relativistic equilibrium, m is the par- 
ticle rest mass, and tm represents a characteristic time 
between subsequent collisions. This can be regarded as 
the relaxation time only in a local rest frame where the 
momentum of the particles is zero[23j]. It is well-known 
that in a general inertial frame, the relaxation time tm 
can be written as follows: 



P 



(10) 



Although, in the Marie model, the transport coefficients 
are expressed usually as functions of the characteristic 
time tm, they cannot be described as functions of the 
relaxation time tm because it depends on the micro- 
scopic momentum component p°, which means on mi- 
croscopic jy— , 1 .) , and therefore it cannot appear 

in any macroscopic description. To avoid this problem, 
Takamoto and Inutsuka [26( proposed a modified Marie 



4 



model, in which the relaxation time r is taken as the 
weighted average, i.e. ^=(^-). With this approxima- 
tion, the following relation can be obtained [IM [2(| 



t~M 



Ki(x) 
K 2 (x) 



(11) 



where x=-^y- and K n is the second kind modified Bessel 
function of order n. The correction ^ 1 i x \ tends to 1 at 
low temperatures, i.e. ?>oo, and to ^ in the limit of 
high temperatures, i.e. >0. In this modified approach, 
the characteristic time in the transport coefficient can be 
replaced by the relaxation time r using Eq. as an 
approximation. 

The Marie model provides a good approximation of the 
full collision term at low temperatures. 

A more general RBGK model, which provides a reason- 
able approximation of the transport coefficients at both 
low and high temperatures, was subsequently proposed 
by Anderson and Witting [23], and it reads as follows: 



C 2 TA 



(f eq - f) 



(12) 



ta being the relaxation time. 

Both models can reproduce on the macroscopic level 
the conservation equations given by d^T^ — 0, and 
d^N^ — 0, although with different expressions for the 
dissipativc terms and transport coefficients. For in- 
stance, the shear viscosity using the Marie model is given 



by J]M r - 



,4P e «7 



for high temperatures (ultra-relativistic 



case), with P eq the equilibrium pressure, while with the 
Anderson- Witting model yields tia— 4P 5 ta ■ 

In general, the dissipation parameters, like the bulk 
viscosity, thermal conductivity and shear viscosity, are 
only approximations of the values obtained by lineariza- 
tion of the full collision term in the relativistic Boltzmann 
equation, EqJT] 

Having discussed the BGK formulation in the relativis- 
tic context, we next proceed to map it within the Lattice 
Boltzmann framework. 



B. Lattice Boltzmann Model 

The Lattice Boltzmann theory for classical fluids shows 
that it may prove more convenient to solve fluid problems 
by numerically integrating the underlying kinetic equa- 
tion rather than the macroscopic fluid equation them- 
selves. The main condition for this to happen is that a 
sufficiently economic representation of the velocity de- 
grees of freedom be available. Following upon con- 
solidated experience with non-relativistic fluids, such a 
representation is indeed provided by discrete lattices, 
whereby the particle velocity (momentum) is constrained 
to a handful of constant discrete velocities, with sufficient 
symmetry to secure isotropy and the fundamental con- 
servations of fluid flows, namely mass-momentum-energy 



conservation and rotational invariance. The main ad- 
vantages of the kinetic representation of classical fluids 
have been discussed at length 28], and they amount basi- 
cally to the fact that the information is transported along 
straight-streamlines (the discrete velocities are constant 
in space and time) rather than along space-time depen- 
dent trajectories generated by the flow itself, as it is case 
for hydrodynamic equations. Moreover, diffusive trans- 
port is not described by second-order spatial derivatives, 
but rather emerges as a collective property from the adia- 
batic relaxation of the momentum flux tensor to its local 
equilibrium value. This is crucial in securing a balance 
between first-order derivatives in both space and time, 
which is essential for relativistic equations. 

In order to reproduce the relativistic hydrodynamic 
equations, an LB model with the D3QI9 (19 speeds in 3 
spatial dimensions) cell configuration, as shown in Fig.[T] 
was proposed in Ref. From Fig. [TJit is readily appre- 



' Si 



be- 



ciated that the highest D3QI9 speed is y/2ci, cr- 
ing the limiting lattice speed along each direction. The 
velocity units are rescaled such that the speed of light 
becomes c=l. 

As noted above, relativistic hydrodynamics evolves two 
scalars, number and energy density. It is therefore con- 
venient to introduce two separate distribution functions 
fi and <7i for each velocity vector c\, representing, so to 
say, "fluons" and "phonons", respectively. 

The hydrodynamic variables are calculated by using 
the following five macroscopic constraints, 



18 



= £/< 



18 



i=0 



18 

(e + P) 7 2 u = ^ 9iCi 

8=0 



(13a) 



(13b) 



(13c) 



From these equations, we need to extract six physical 
quantities, namely n, u, e and P. With five equations 
for six unknowns, the problem is closed by choosin g an 
equation of state, which we take of the form e=3P [23| . 
We wish to emphasize that the present LB scheme is by 
no means limited to this choice. 

Both distribution functions fi and gi are postulated 
to evolve according to the relativistic Boltzmann-BGK 
equation based on the low-temperature Marie model, 
Eq. ©. 

To obtain the lattice analogue of the Marie model, we 
first write explicitly Eq. © as follows: 



do(p°f)+d a (p a f) 



tm 



-(/ eq -/) 



(14) 



Replacing the value of the four-momentum, we obtain 



do(mj v f) + d a {mc a -i v f) 



TM 



-ir-f) 



(15) 
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with 7„ the Lorentz factor for the microscopic velocities 
c a . Due to the fact that the velocity and spatial coordi- 
nates are linearly independent, we can further write: 



Jvdof + JvC a d a f 



f cq -f 



(16) 



Dividing by j v on both sides of Eq. (|16l) , we obtain 



dof + c a d a f 



tmP 



o(n-f) 



and replacing Eq. (fl"0|) . we obtain 

1 



dof + c a 8 a f 



(f eq - /) 



(17) 



(18) 



According to the modified Marie model [26j , we can write 
Eq. GHJ) as 



dof + c a d a f 



<? (/ cq -/), 



(19) 



where the correction term using Eq. (|lip . is given by 

-)"- 

tm* I tm* 
'Kx{ X ) 1 



1 



(20) 



For low temperatures, ^[^j 



1 and 7„ ~ 1, so that the 



correction term $ tends to zero, thereby renstituting the 
non-relativistic Boltzmann equation. At high tempera- 
tures, this term can be approximated by 



a-* J- 

2 T M 



(21) 



which also tends to vanish as temperature is made higher. 

As noted above, Eq. (fT9|) . without the term is just 
the Boltzmann equation for the case of non-relativistic 
fluids, with the collision time r representing a realistic 
relaxation time of the system. 

Therefore, for the purpose of this work, we postulate 
the discrete distribution functions to evolve according to 
the following pair of BGK Boltzmann equations p9| . 

fi (x + CiSt,t + St)-fi (x, t) = - f (fi - jf 1 ) , (22) 

and, 

9i (x + CiSt,t + St)- 9l (x , t) = - f ( 9i - g?) , (23) 

where f° q and g eq are the equilibrium distribution func- 
tions. 

To find the equilibrium distribution functions recover- 
ing the relativistic fluid equations, Eqs. ([5]) and ©, in 
the continuum limit, we use the moment-matching pro- 
cedure described earlier on in this paper. 

More precisely, we write the equilibrium distribution 
functions as, 



fP = Wi [A + Ci-B] ,for i>0 



(24a) 




FIG. 1: Set of discrete velocities for the relativistic lattice 
Boltzmann model. The highest speed is V2ci. 



g ° q = Wi [C + Ci ■ D + E : (cicl - al)} Tor i>0 



gl q = Wo[F} , 



(24b) 
(24c) 



where a, A, B, C, D, E and F are Lagrange parameters, 
to be fixed by matching the discrete to the correct con- 
tinuum equations. The weights Wi for this set of discrete 
speeds are defined by wo = 1/3 for the rest particles, 
Wi = 1/18 for the velocities |ci|=Q, and Wi = 1/36 for 
\ci\ = V2ci. 

First, we find the values for A and B to obtain the 
conservation of particle number, Eq. ©. To this purpose, 
we impose 



18 



(25) 



2 = 



and, 



in 



i=0 



■en -> -> 

. H a = n-fu 



(26) 



Replacing the Eq. (f2"4"]l into the sums, we arrive to 



18 



and 



18 



(27) 



(28) 



where, we can see easily that A=wy and B=\'ynu. Next, 

we have to obtain the Eq. (JS|) from the equilibrium dis- 
tribution functions g cq . To this end, we impose the fol- 
lowing constraints: 



18 



(29) 
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18 

5>r^ = (e + P) 7 2 u 

i=0 



(30) 



and additionally, 

18 



^ fJ^CiaCip = PS cl b + (e + P)j 2 U a U b 



(31) 



i=0 



Using a similar procedure as before, we can find 
the rest of the Lagrange parameters, ct=-g-, C , = ^f, 
D=\{e + P) 7 2 u, E ab ^(e + P)j 2 u a u b , and P=(e + 



P) 7 2 



3-3 



(2+<)P 



^(e + P) 7 2 |up 



These calcu- 



lations are shown in detail in Appendix [Al 

The equilibrium distribution functions recovering the 
relativistic fluid equations in the continuum limit, finally 
read as follows: 



eq _ 



WiWy 



1 



, (3 • u) 



(32) 



for i>0, 



^ = w t (e + Ph 2 



3P 



, (Cj • u) 



(P + e) 7 2 c 2 
9 (gj u) 2 



3H 2 
2 c 2 



(33) 



for i>0, and 



9o^ 



w {e + P)j 2 



3- 



3P(2 + c 2 ) 3|u|' 



(P + e) 7 2 c 2 



2 c 2 



(34) 



for the rest particles. 

By Taylor expanding the Eqs. (f2"2"]l and (f2"3")l to second 
order in <5t, and retaining terms only up to first order 
in the Chapman-Enskog expansion / = f eq + nf 1 + . . . , 
where k ~ crV is the Knudsen number, the LB equations 
can be shown to reproduce the following continuum fluid 
equations as derived in detail in Appendix [Bl 



ft [(C + P) 7 2 - P] + 9a [(C + P)7 2 "a] = , 

ft [(c + P) 7 V] + ftP + ft [(e + P) 7 2 U a Ub 
= ft [ft(?77 u a) + d a (r]ju b ) + di{ri^ui)5a b } , 

for the energy momentum conservation, and 
ft(n 7 ) + ft (n-yua) = , 



(35a) 



(35b) 



(36) 



for the conservation of particle number. The indices a,b 
and I denote the spatial components. 

The choice of the state equation, e=3P, simplifies the 
equilibrium functions as follows, 



n 



1 



(37) 




BAMPS: n/s =0.001 — 
BAMPS: n/s =0.05 
BAMPS: n/s =0.1 — 
LBS: n/s =0.001 — 
LBS: n/s =0.05 - 
LBS: n/s =0.1 — 
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FIG. 2: Comparison between the BAMPS simulations [211] and 
the lattice Boltzmann results at £=3.2fm/c. Note that, in 
both simulations, the value of /3~0.2 for the speed of prop- 
agation of the shock wave is obtained. Pressure (top) and 
velocity (bottom) of the fluid as function of the spatial coor- 
dinate z. 




for i>0 and, 



m;e 7 



for i>0 and, 



(cj • u) (&i ■ u) 2 



7 2 c 2 



6- 



eq 2 

9a = wo £7 



4- 



7 2 c 2 



- V- 



(38) 



(39) 



for i=Q. Then, the equations for the macroscopic vari- 
ables take the form: 



4 2 -> 



"7 = fe(7 2 -|) = 

^ i=0 (7iCi. The shear viscosity 



Eifosf and 

is computed as ?7=| 7 e(T — 6t/2)cf. 

Also, it is worth noting that our scheme smoothly re- 
covers the non-relativistic limit by simply letting /3 — > 0. 



IV. DISSIPATIVE HYDRODYNAMICS 

According to kinetic theory, dissipative effects emerge 
at the level of first order terms in the Knudsen number 
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Initial condition - 
BAMPS: t = 0.64 fm/c ■ 
BAMPS: t = 3.20 fm/c 
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FIG. 3: Time evolution of the shock wave for BAMPS 
simulations [2ll] and Lattice Boltzmann results. Here the 
speed of propagation of the shock wave /3 ~ 0.2 is obtained. 



G-e Grid points: 50 
b-h Grid points: 400 
^>Grid points: 800 
^4 Grid points: 3200 
— BAMPS 




FIG. 4: Lattice Boltzmann simulation of the shock wave for 
different grid resolutions, /3~0.2 and 77/s=0.05. 



expansion of the kinetic equations. At a more fundamen- 
tal level, dissipation is an emergent property resulting 
from the finite-time relaxation of non-equilibrium kinetic 
excitations on top of the hydrodynamic "ground state". 
A detailed Chapman-Enskog analysis (see Appendix [Bj , 
shows that the lattice formulation needs to retain second 
order terms in the lattice spacing, which means that the 
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FIG. 5: Relative convergence error for different grid res- 
olutions as a function of the ^-coordinate for /3~0.2 and 
r)/s=0M. 



n n is = o.oi 
O n is = 0.05 

slope = -1.92 
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Log 1Q (A/50) 



FIG. 6: Relative convergence error as a function of the num- 
ber of grid points, for /3~0.2. A represents the number of 
grid points. Here, the relative error E ave is calculated by 
taken the mean value of the relative errors at every location 
z E -AV A 



=1 E r (z). 



streaming operator needs to be expanded to second order 
in the lattice time step 5t and, by the light-cone rule, in 
Sx too. Straightforward but lengthy algebra, leads to the 
following expression of the LB dynamic viscosity 30] 



V = PC S 



5t 
2 



9<?s T 



1 - - 
2 



where we have defined 



St 



(40) 



(41) 
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as the parameter measuring the time-granularity of the 
LB scheme. Indeed, the limit 9 — > reproduces the con- 
tinuum value r\ = pc 2 s T. Similar calculations for the rela- 
tivistic case yield (see more details in Appendix [Bj 



77 = 7cMe + P)(l-0/2)c?. 



(42) 



A few comments are in order. First, we note that pos- 
itivity of the kinematic viscosity implies 



< 9 < 2 



(43) 



This linear stability constraint for the discrete scheme, is 
readily seen to associate with the second-principle (neg- 
ative viscosity implies physical instability). 

The above expressions seem to suggest that ideal hy- 
drodynamics, i.e. strictly zero dissipation, could be 
achieved in the limit St —> 2r, i.e. 9^2. Actual 
practice, though, shows that this limit is an illusory 
one, since, whenever the viscosity falls below a given 
(flow-dependent) threshold, the stability of the scheme 
is compromised. Physically, the reason is that below a 
given threshold, the system is no longer capable of dissi- 
pating short-scale gradients, thereby allowing the non- 
equilibrium component of the distribution function to 
grow wildly, and finally ruin the simulation. This is 
in line with the so-called "numerical uncertainty prin- 
ciple" (NUP) for transport advection equations, accord- 
ing to which a minimum non-zero viscosity is required to 
secure the positivity of the positive-definite quantities, 
such as the fluid density 31|. In a nutshell, the point is 



that, in order to reach zero viscosity with a positive def- 
inite distribution, wavelengths at all scales are needed, 
including those below the lattice spacing Sx. Since - 
by construction- the latter are missing from a discrete 
lattice representation, positivity can only be maintained 
through a finite amount of dissipation, typically of the 
order of the inherent lattice viscosity . Incidentally, 
we note that viscosity has the same physical dimension 
as SxSv ~ h/m, whence the notion of "uncertainty prin- 
ciple". 

For LB equations, the NUP can be formulated in 
terms of an inequality involving the equilibrium and non- 
equilibrium components of the discrete distribution func- 
tion. To appreciate this point, let us first recast the stan- 
dard LB in the following collide-stream form: 

Mx + c^t + St) = f[{x-t) = (l-9)f l (x;t) + 9f^ (44) 

where /' denotes the so-called post-collisional distribu- 
tion function. 

From the above, it is seen that positivity of the post- 
collisional distribution at time t guarantees positivity of 
the distribution at the subsequent time t + St. Simple 
algebra yields: 

l/fl 



9 < 9 NUP [f] = min 



If, 



neq i 



1. Weak non-equilibrium (9nup > 2) 

2. Strong non-equilibrium (1 < 9nup < 2) 



3. Extreme non-equilibrium 



7 NU P 



<1) 



This informative expression suggests the definition of 
three distinct non-equilibrium regimes: 



In the weak non-equilibrium regime (often referred to 
as strong-coupling regime) , the one relevant to hydrody- 
namics, the NUP does not set any additional constraint 
to linear stability. In the strong non-equilibrium regime, 
however, non-linear stability may in principle set the 
most stringent constraint. Clearly, this is even more so 
in the extreme non-equilibrium region, where the non- 
equilibrium component exceeds the equilibrium one, in 
total defiance of hydrodynamics. 

Remarkably, LB proves capable of stable operation in 
this "linearly-forbidden" region. In fact, the negative 
shift —St/2, ( "propagation viscosity" in LB jargon) which 
stems directly from the light-cone structure of the LB 
streaming operator, permits to attain very small viscosi- 
ties, of order, say, 10~ 3 in lattice units, while still keeping 
St = 1, and 9 ~ 2 - 0(1CT 3 ). This allows for the simu- 
lation of very-low viscous flows (such as the quark-gluon 
plasma) with time-steps of order O(l), which proves very 
beneficial for computational purposes. 

The ultimate reason for such favorable behavior in the 
strong non-equilibrium regime can be traced to the exis- 
tence of lattice versions of the H-theorem [32|, HU . 

Another remarkable property of the LB formulation 
is that, in contrast to hydrodynamic formulations, dis- 
sipation is not represented explicitly through second- 
order spatial derivatives, but emerges instead from a 
first-order, covariant propagation-relaxation dynamics, 
through adiabatic enslaving of the momentum- flux tensor 
to its equilibrium (ideal-hydrodynamic) expression. As a 
result of this first-order dynamics, the CFL (Courant- 
Friedrichs-Lewy) stability condition of the LB scheme 
reads simply as uSt < Sx, instead of vSt < Sx 2 , the latter 
being much more demanding on the time-step St, as the 
grid is refined (Sx — > 0). In the above, v — is the fluid 
kinematic viscosity. 

Also to be noted, built-in causality is secured by the 
hyperbolic structure of the underlying kinetic theory. 

Before closing this section, we wish to emphasize that 
the structure of the dissipative terms could be enriched 
by turning to a multi-time relaxation version of the colli- 
sion operator, whereby different moments relax with dif- 
ferent rates to their equilibrium expression [34], 35] . This 
allows to enlarge the list of transport coefficients, includ- 
ing bulk viscosity, thermal conductivity and anisotropic 
transport parameters. 



V. VALIDATION AND APPLICATIONS 

Having discussed the basic aspects of the relativistic 
Lattice Boltzmann theory, we next move on to its numeri- 
cal validation and application to two different problems of 
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modern relativistic hydrodynamics, namely shock prop- 
agation in viscous quark-gluon plasmas and blast-waves 
from supernova explosions in interstellar media. 



A. Quark-Gluon Plasma 

To test the model, we solve the Riemann problem in 
viscous gluon matter [2l| with a ultra- relativistic equation 
of state e=3-P, as before, and the relation between energy 
density and par ticle number density, e=3nT, T being the 
temperature [23J. The initial configuration consists of two 
regions, divided by a membrane located at z=0. Both 
regions are thermodynamically equilibrated, at different 
constant pressure, Pq for z<0 and Pi for z>0. At t=0, 
the membrane is removed and the fluid starts expanding. 

We implement a one-dimensional simulation with an 
array of size 1 x 1 x 800 using open boundary conditions 
at the two ends of this ID chain. In this case, the 4- 
velocity is given by u f± =(j, 0, 0, 7/3) M . The velocity of the 
lattice is chosen q = 1.0, so that the cell size Sx and time 
step St are both fixed to unity. This corresponds in IS 
units to (5a;=0.008fm and i5t=0.008fm/c. The viscosity is 
calculated as T}=^e(r — 1/2), and the entropy density by 
the approximation s=4n — nlnA, with A=^Sy the gluon 
fugacity. The equilibrium particle density n eq is given 
by, n eq = dc ^ with c?g=16 for gluons. Next, we calculate 
the ratio between the viscosity and entropy density, rj/s, 
that is used as a parameter to characterize the conditions 
for the onset of shock- waves. The pressures were chosen 
as P =5.43GeVfm- 3 and Pi=2.22GeVfm" 3 , correspond- 
ing to 7.9433 xl0~ 6 and 3.2567xl0~ 6 in numerical units, 
respectively. The initial temperature is Xb=350MeV, cor- 
responding to To=0.0287 in numerical units. With these 
parameters, the conversion between physical and numer- 
ical units for the energy, is lMeV=8.2x 10~ 5 . 

Fig. [2] shows the results for different values of r\j s and 
the comparison with the B AMPS [36] (Boltzmann Ap- 
proach of Multiparton Scattering) microscopic transport 
model simulations [2lj at time 3.2fm/c. Fig. [3l shows 
the time evolution of the system for ?7/s=0.1 for the two 
numerical models. In both cases, excellent agreement 
with BAMPS is observed. Fluids moving at higher speed, 
/3~0.6, were also considered in Ref. [l|, where numerical 
"tachyons" with c;=10 were used. 

Indeed, from Eqs. (f3"2"j) and (j3"3")l . we see that the pos- 

itivity condition, /? q >0, implies Cj • u<-§-. As a result, 
by raising cj, e.g. by reducing the time-step accordingly, 
positivity can be preserved for higher values of (3. 

To check the convergence of the model, we implement 
simulations taking rj/s=0.05 and 77/3=0. 01 for different 
grid resolutions. Fig. |4] reports the pressure profile at 
time 3.2fm/c and shows very small differences between 
the results when the resolution is changed from 50 to 
3200 grid points with 77/s=0.05. To obtain a more quan- 
titative measure of the convergence we use the Richard- 
son extrapolation method [13, In this method, given 
any quantity A(Sx) that depends on a size step Sx, we 



Grid points 


Total time steps 


CPU time (ms) 


50 


25 


0.94 


100 


50 


3.5 


200 


100 


17.1 


400 


200 


68.4 


800 


400 


272 


1600 


800 


1095 


3200 


1600 


4396 



TABLE I: Computational time required for the simulation of 
the shock waves in quark-gluon plasma as a function of the 
grid resolution. 



can make an estimation of order n of the exact solution 
A by using 



A = lim A(5x) 

Sx— ¥0 



2 n A (4f ) - A(Sx) 



+ 0(Sx n+L ) 



(45) 

with errors 0{Sx n+1 ) of order n + 1. Thus the relative 
error between the value A(5x) and the "exact" solution 
A can be calculated by 



E r {5x) 



A{Sx) - A 



A 



(46) 



In our case, the quantity A is the pressure P(z) and we 
set up n = 2. We can estimate the relative error as shown 
in Fig.[5]for ?7/s=0.01 using Eqs. (|45l and (|46|) . at every 
grid point. Indeed, the relative error with respect to the 
"exact solution" decreases rapidly with increasing grid 
resolution. More precisely, Fig. [5] shows that the present 
scheme exhibits a near second-order convergence. This is 
basically in line with the convergence properties of non- 
relativistic LB schemes. 

However, we can see that for higher viscosity, i.e. larger 
values of the relaxation time r, and higher grid reso- 
lution (smaller Sx), the order of convergence decreases 
due to the lack of adiabaticity associated with increasing 
Knudsen number. Nevertheless, the model is still able to 
reproduce shock waves, at low resolution, hence with a 
very modest computational time. For instance, using a 
resolution of 50 grid points, the simulation took 0.94ms 
in a standard PC. Other values are shown in Table |U 
From this table, it is readily appreciated that the com- 
putational cost scales linearly with the number of grid 
points and time-steps. 



B. Supernova explosion simulation 

Several important astrophysical phenomena involve 
strongly-relativistic hydrodynamics, and some of them 
fall in the region of 7 ~ 1.4, covered by our scheme. This 
is the case, for instance, of blastwaves produced by super- 
nova explosions [20] . In this section, we simulate a shock 
wave, generated by, say, a GRB (7-ray burst) or XRF (X- 
ray flash) supernova explosion [19i . |20| , colliding against 
an interstellar cloud composed by massive matter, e.g. 
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(a) 




(b) 




(c) 



FIG. 7: Relativistic shock wave, generated by a 7-ray burst 
or X-ray flash supernova explosion [l|| [2(j, impacting on a 
massive interstellar cloud at \/3\ — 0.5 at t = 1350 time steps 
which is equivalent to 4280 years. Here the streamlines repre- 
sent the velocity field, and the colors (a) the pressure, (b) the 
particles density, and (c) the temperature. The simulation 



was implemented on a grid of 200 x 100 x 100 cells. 



molecular gas[17]. The ejecta from the explosion of such 
supernovae are known to sweep the interstellar material 
along, up to relativistic velocities (relativistic outflows) 

The simulation is implemented in a box of size 6x3x3 
xlO 16 Km in a coordinate system (x,y,z), using a lat- 
tice of 200 x 100 x 100 cells, which gives a cell length 
5x=5y=5z=3 x 10 14 Km, using numerical "tachyons" 
with c;=10, and a time step St=3.17 years. The simula- 
tion region is divided in two zones by the plane x = 50. 
The interstellar medium, located at x > 50, is character- 



ized by a particle density ni=0.6 cm -3 and temperature 
Ti=10 4 K. The massive cloud is modeled as a spheri- 
cal obstacle, with a radius of 10 cells, centered at loca- 
tion (100,50,50). The boundary condition on the sur- 
face of the obstacle is implemented forcing the obstacle 
cells to evolve to the equilibrium distribution function 
with the constant values, n = n\, u = 0, and T = T\. 
Open boundary condition was implemented at right, left, 
top, bottom and front of the simulation zone according 
to the shock wave propagation direction (x-direction) , 
which consists on copying the information of the distri- 
bution functions from the second last cells to the last 
ones of the boundary. At back boundary we set an in- 
let flow boundary condition fixing the distribution func- 
tions of the boundary cells with the equilibrium distri- 
bution function evaluated with the initial conditions no 
and To (30L l39j . In order to obtain a shock wave moving 
at \f3\ ~ 0.5 along the x-direction, we set Tq=QT\ and 
riQ=1n\ for the region x < 50. The simulation, span- 
ning 1350 time steps, takes about 1900 CPU seconds on 
a standard PC. Fig. [7] shows the simulation results for 
the velocity, pressure, particle density, and temperature 
fields of the supernova remnant, during the impact of 
the shock wave on the massive interstellar cloud, red and 
blue denoting high and low values, respectively. 

Here, we can see that the density n is higher in the 
shock front, due to sweeping of interstellar material by 
the shock-wave, which is compressing the fluid. On the 
other hand, the temperature of the fluid is higher in the 
zone of x < 50, as a consequence of the initial configu- 
ration. The temperature is seen to increase in the zone 
where the collision takes place (see Fig. |5J), and so does 
the temperature. This is due to conversion of kinetic en- 
ergy to pressure/ temperature caused by the momentum 
lost on the solid boundary of the massive cloud. 

Fig. [8] illustrates in more detail the density n, pres- 
sure P and temperature T of the fluid during the col- 
lision and compares the respective curves with the ones 
obtained when the obstacle is absent. Note that the par- 
ticle density, pressure, and temperature values, with and 
without obstacle, present a small difference sufficiently 
downstream the obstacle along the cc-axis at y — z — 50 
(see Fig. [8|. During the collision, the Shockwave sur- 
rounds the obstacle and later the fluid meets again at the 
x-axis and overlaps. Due to this, the x-component of the 
Shockwave propagation velocities are the same (because 
of symmetry) for all the incoming fluid to the meeting 
zone, the perturbations along this axis close to the shock 
front are weak, contrary to the zone near the obstacle, 
where the fluid fills up again, due to the low pressure. 
However, the fluid moves slower than in the case without 
obstacle because of the existence of flow moving outwards 
off the axis. If we increase the ratio between the cross 
section and the length of the obstacle, larger departures 
between the velocity of the shock-fronts with and without 
obstacle would be expected. Moreover, later in time after 
the collision, differences in the pressure and other quanti- 
ties, can generate turbulence. Transversal perturbations 
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FIG. 8: Pressure P, number of particles density n, and tem- 
perature T of the supernova remnant as a function of the x 
coordinate at y = z = 50 and t = 1350 time steps equivalent 
to 4280 years. The solid lines describe the values in the pres- 
ence of the obstacle, dotted line the region where the massive 
interstellar cloud is located, and the dashed line the values 
without obstacle. 



in the variables, as one moves out from the x-axis, are 
shown in Figs. M and HUJ 

Shock waves form when the speed of injection of mass 
exceeds the sound speed of the surrounding medium [l]} ■ 
By changing the values of the temperature of the fluid in 
the region x < 50, in order to obtain speeds of mass in- 
jection of |/3| = 0.5, = 0.2, and |/3| = 0.01, we can 
see that the increment of the particle density due to the 
sweeping of interstellar medium by the shock wave be- 
comes appreciable only for |/3| = 0.5 (see Fig.QTJ]). A sim - 
ilar argument applies to the pressure cone (see Fig.[ ^a)| . 
Indeed, in the other cases, the speed of mass injection is 
lower than the sound speed, and therefore no shock-wave 
can be formed. 
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FIG. 9: Fluid pressure after the collision of the shock wave, 
produced by a supernova explosion, against the massive inter- 
stellar cloud, at [(a)] | /3 1 = 0.5, [(b)] |/3| = 0.2, andO|^| = 0.01 



in a cut going through the center of the cloud. The stream- 
lines represent the velocity field, and the colors the pressure. 
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FIG. 10: Particles density of the fluid after the collision of 
the shock wave, produced by a supernova explosion, against 



the massive interstellar cloud, at (a) \/3\ = 0.5, (b) = 0.2, 



and (c) j/3| = 0.01 in a cut going through the center of the 
cloud. The streamlines represent the velocity field, and the 
colors the particles density. 



VI. CONCLUSIONS AND OUTLOOK 

In this paper, we have provided a detailed discussion of 
the Lattice Boltzmann formulation for relativistic fluids. 
In particular, details on the construction of the relevant 
lattice equilibria are provided, emphasizing the common 
aspects with standard Lattice Boltzmann theory. 

The scheme is shown to exhibit excellent agreement 
with previous numerical simulations of shock wave prop- 
agation in quark-gluon plasmas, at a fraction of the cost 
of hydrodynamic codes. Near-second order accuracy with 
grid resolution and linear computational time with space- 
time resolution, are evidenced. 

As an example of relativistic hydrodynamics with non- 
trivial geometries, we have also applied our scheme to 
an astrophysical system, namely the collision of a shock 
wave, produced by a supernova explosion, against a cold 
molecular cloud. The numerical simulations show good 
qualitative results yielding information, that can be com- 
pared with experimental results and other numerical 
methods. 

For the case of quark-gluon plasma simulations, the 
present lattice-kinetic algorithm appears to be nearly an 
order of magnitude faster than corresponding hydrody- 
namic codes. This is due to the fact that, at variance 
with any hydrodynamic representation, LB moves infor- 
mation along constant light-cones rather than space-time 
changing material fluid streamlines (ifjj . This trivializes 
the Riemann problem to a mere shift of the distribution 
function along the corresponding lightcone, a floating- 
point free, exact operation, which is way more convenient 
than propagating hydrodynamic fields along space-time 
changing streamlines. Such an advantage, key in ordi- 
nary lattice Boltzmann fluids, might be even accrued in 
the relativistic context. 

Several issues remain open for future research. First, 
extensions of the present scheme to higher-order lattices 
are worth being considered, for they should give access 
to higher values of /3, by use of correspondingly higher- 
order lattice equilibria. This strategy has indeed proved 
very effective for the case of compressible and thermal 
non- relativistic fluids [4ll - l44| . 

Another important question concerns the existence of 
a relativistic lattice H-theorem. Apart from the theoret- 
ical interest on its own, this has major implications on 
the numerical stability of the scheme at high Reynolds 
number, i.e. for the simulation of relativistic turbulence 

m. 

Yet another interesting research direction is the simu- 
lation of relativistic flows with a non-ideal equation of 
state, which may find applications in relativistic cos- 
mology and high-energy theories of the early universe 

EE S3- 

These are just but a few of the many exciting develop- 
ments and applications which may currently be envisaged 
for the relativistic Lattice Boltzmann equation presented 
in this paper. 
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Appendix A: Moment Matching Procedure 

To obtain the equilibrium distribution functions f^ q 
and g^ q that reproduce in the continuum limit the hy- 
drodynamic equations, Eqs. (J5J) and ©, we use the 
moment-matching procedure. In section MI B[ we de- 
scribe the procedure and calculate the equilibrium distri- 
bution functions /° q in order to obtain the conservation 
of particle number, Eq. ([5]). Following a similar proce- 
dure, to find the equilibrium distributions g° q , first we 
can write it, as before, as 



g^ = Wi[C + Ci-D + E: (cici-al)] ,for i>0 



g c Q q = w [F] 



(Ala) 
(Alb) 



with C, D, a, and E the Lagrange multipliers. Then, we 
impose the following constraints: 



is 



is 



and additionally, 
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gt q CiaCi(3 = P8ab + (e + P)^ 2 U a U b 



(A2) 



(A3) 



(A4) 



Replacing Eq. ([Ml) into Eq. ([Ml) . ([Ml) , and (|A4|) . and 

summing up over the index i, we obtain 

i(2C + F + Tr(2)(c?-2a)) = 1 2 (e + P)-P , (A5) 



and 



C j-D=(e + P) 1 2 u , 



2 

-i- [3C + (cf - 3a)Tr(2)) 5 ab + -^E ab = PS ab 

+ (e + P)j 2 u a u b , 



(A6) 



(A7) 



where we have defined Tr(E) as the trace of the tensor E. 
From Eq. (|A6|) we can see that D=^(e + P)^ 2 u. If we 



compare the left and right hand sides of Eq. (|A7[) . we can 

2 

conclude that a=^-, and therefore Eq. (|A7[) is simplified 



to 



-CS ab 



2r 4 



PSab 



P)l u a u b 



(A8) 



Comparing again both sides of this equation the La- 
grange multipliers C=^ and E ab =^(e + P)^ 2 u a u b are 
obtained. Now, the only missing parameter to be deter- 
mined is F. Replacing the values of C, a, and E into 
Eq. ||A5]) . it gives 



2-P F cf o 

~cj + 1 + i Tr{E) = 7 (e + p) 



(A9) 



From here, we can get the Lagrange parameter F and it 
can be written as 



F = (e+P)Y 



3-3 



(2 + cf)P 
c 2 {e + P)j 2 



^ + Ph 2 \u\ 



(A10) 

Summarizing, we have determined all the Lagrange pa- 
rameters and therefore the equilibrium distribution func- 
tions g eq that recover in the continuum limit the conser- 
vation equation for the momentum-energy. 



Appendix B: Chapman-Enskog Expansion 

The discrete Boltzmann equations, Eqs. (f2"2"j) and (|23l) . 
determine the evolution of the lattice relativistic fluid. 
In the continuum limit, these evolution rules must re- 
produce the partial differential equations of relativistic 
hydrodynamics. In order to accomplish this task, we 
adopt a standard Chapman-Enskog expansion. We start 
by taking the Taylor expansion of the Boltzmann equa- 
tions, up to second order in spatial and temporal coordi- 
nates, 



1 



a.b 

+ dtvMi + \d 2 lM 2 = -~(fi - /D 

l T 



Viad a gi + d ad b giVi a v ib + d t g l 

a,b 

+ d t v ia d a gi + \d 2 g t 8t 2 = --{g l - g° q ) , 

Z T 



(Bla) 



(Bib) 



where a, b=x, y, z denote the x, y and z components. 
Next, we expand the distribution functions, and the 
space-time derivatives in a power series of the Knudsen 
number follows: 



(B2a) 
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(0) , (1) . 2 (2) . 

9i = 9\ + K 9i + K 9i + 



d t = Kd tl + n 2 d t . 



(B2b) and 
(B2c) 



d t2 (nj) + d 2a (nju a ) = 



(BIO) 



da = Kdla + K 2 d2a- 



(B2d) 



It is assumed that only the Oth order terms of the distri- 
bution functions contribute to the macroscopic conserved 
variables. Therefore, for n > we have 



E^ = ° > E^ = ° • 



(B3a) 



(B3b) 



By inserting these results into Eqs. (|Blal) and (IBlbl) . 
we obtain at Oth-order in n 



oq (0) 

9i =91 



to the hrst order in k, 



p(0) 



f(l) 



T 



(B4) 



(B5) 



(B6a) 



(B6b) 



and to the second order in k, 



d t j o) +v ia d 2a fr = - 



(0) 



(2) 



(B7a) 



- ^ (W 4 a5la +9 fl )ff- 1) 

(2) 

+ Ot 2 9i +Via02a9 l - 



(B7b) 



A this stage, all the ingredients required to deter- 
mine the equations that the model satisfies in the con- 
tinuum limit, are available. By summing up Eqs. (|B6a|) . 
(IB6bl) . (|B7ap . and (|B7b|) over index i, taking into ac- 
count Eqs. (|B4[) . (|B5I) . and the equilibrium distribution 
functions defined by Eqs. (|3"T)l . (1551) . and (|39p . we obtain 



d tl (nj) + diaiwyua) = , 



d tl ((e + Ph 2 -P) + d la ((e + P) 1 2 u a ) =0 , (B9) 



d t2 ((e + P) 1 2 -P) + d 2a ((e + P) 7 2 u a ) = . (Bll) 

By adding these equations, the first and second scalar 
equations, associated with the conservation of the num- 
ber of particle and the first conservation equation for the 
momentum-energy, 



dt{nrf) + da^wyua) = 



(B12) 



and 



d t ((e + Ph 2 - P) + d a {{e + Ph 2 u a ) =0 , (B13) 



are obtained, which correspond to Eqs. (|36l) and (|35a[) . re- 
spectively. To derive the second conservation equation, 
Eq. (|35b[) . the equations (|B6b|) and ()B7b|) must be multi- 
plied by Vi and summed up over the index i, which leads 
to 



+ dia [(e + P)j 2 u a u b ] 



(B14) 



and 



dt 2 [(e + P)-f 2 u b ]+d 2b P 

+ d 2a [(e + Ph 2 u a u b ] + duES! = 



(B15) 



where the first order tensor n 



(i) 

ab 



(i 



(1) 



(1) 



is defined. By replacing the distribution function / 

from Eq. (|B6b[) into the tensor H^ b \ and the result into 
Eq. (|B15|) . we obtain 

d t2 [(e + P) 7 2 u fc ] + d 2b P + 8 2a [(e + Ph 2 u a u b ] 

- di a [dib{viu a ) + d la (riju b ) + du(rrfUi)5 ab ] = , 

(B16) 

with the viscosity 77=^7(6 + P)(t — dt/2)cf, I denoting 
again the spatial components. To arrive to these re- 
sults, we have assumed low-speed, \u\ <C c. The sec- 
ond momentum-energy conservation equation, Eq. (|35bl) . 
is obtained by summing up Eqs. (|B16l) and (|B14|) . It gives 



d t [(e + Ph 2 u b ] + d b P + d a [(e + P)l 2 u a u b ] 
- d a [d b {rj^u a ) + d a (r/7 u b) + di{rj-fui)5 ab ] = 



(B17) 



The derivation of the dissipative term associated with the 
viscosity r), in Eq. (|B17[) . is obtained assuming low values 
(B8) °f P to neglect higher order terms (~ |u| 3 ) contributions. 

Summarizing, Eqs. (IB12j) . (IB13I) and (|B17I) determine 
the evolution of the fluid, according to the relativistic 
hydrodynamics equations. 
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